#include "common.h"
#include "Interp.h"

void interp (const int iaer, const int idatmp, 
			 const float wl, const float taer55, 
			 const float taer55p, const float xmud, 
			 InterpStruct& is)
{
/*     that for the atmosphere :
     the reflectances
                     rayleigh                             = rorayl
                     aerosols                             = roaero
                     mixing                               = romix
     the downward transmittances
                     rayleigh                             = dtotr
                     aerosols                             = dtota
                     total                                = dtott
     the upward transmittances
                     rayleigh                             = utotr
                     aerosols                             = utota
                     total                                = utott
     the spherical albedos
                     rayleigh                             = asray
                     aerosols                             = asaer
                     total                                = astot
     the optical thickness of total atmosphere
                     rayleigh                             = tray
                     aerosols                             = taer
     the optical thickness of the atmosphere above the plane
                     rayleigh                             = is.trayp
                     aerosols                             = taerp
     the tsca of the aerosols (god dammed it)
                     total atmosphere                     = tsca */
      
	int linf = 0;
	for(int i = 0; i < 9; i++) if(wl > sixs_disc.wldis[i] && wl <= sixs_disc.wldis[i+1]) linf = i;
	
    if(wl > sixs_disc.wldis[9]) linf = 8;
	int lsup = linf + 1;


	/*    interpolation in function of wavelength for scattering
	    atmospheric functions from discrete values at sixs_disc.wldis */
 
	float alphaa = 0;
	float betaa = 0;
	float alphar = 0;
	float betar = 0;
	float alphac = 0;
	float betac = 0;
	is.phaa = 0;
	is.roaero = 0;
	is.dtota = 1;
	is.utota = 1;
	is.asaer = 0;
	is.taer = 0;
	is.taerp = 0;
	float coef = (float)log(sixs_disc.wldis[lsup] / sixs_disc.wldis[linf]);
    float wlinf = sixs_disc.wldis[linf];

    if(iaer != 0)
	{
		alphaa = (float)(log(sixs_aer.phase[lsup] / sixs_aer.phase[linf]) / coef);
		betaa = (float)(sixs_aer.phase[linf] / pow(wlinf,alphaa));
		is.phaa = (float)(betaa * pow(wl,alphaa));
	}

	float d2 = 2 + delta;
    is.phar = (2 * (1 - delta) / d2) * .75f * (1 + xmud * xmud) + 3 * delta / d2;
	if(idatmp == 0)
	{
		betar = 0;
        betaa = 0;
        betac = 0;
	}
	else
	{
		if(sixs_disc.roatm[0][linf] < 0.001)
		{
			is.rorayl = sixs_disc.roatm[0][linf] + (sixs_disc.roatm[0][lsup] - sixs_disc.roatm[0][linf])
				   * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
		}
		else
		{
			alphar = (float)(log(sixs_disc.roatm[0][lsup] / sixs_disc.roatm[0][linf]) / coef);
			betar = (float)(sixs_disc.roatm[0][linf] / pow(wlinf,alphar));
			is.rorayl = (float)(betar * pow(wl,alphar));
		}

		if(sixs_disc.roatm[1][linf] < 0.001)
		{
			is.romix = sixs_disc.roatm[1][linf] + (sixs_disc.roatm[1][lsup] - sixs_disc.roatm[1][linf])
				  * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
		}
        else
		{
			alphac = (float)(log(sixs_disc.roatm[1][lsup] / sixs_disc.roatm[1][linf]) / coef);
			betac = (float)(sixs_disc.roatm[1][linf] / pow(wlinf,alphac));
			is.romix = (float)(betac * pow(wl,alphac));
		}

		if(iaer != 0)
		{

			if(sixs_disc.roatm[2][linf] < 0.001)
			{
				is.roaero = sixs_disc.roatm[2][linf]+(sixs_disc.roatm[2][lsup] - sixs_disc.roatm[2][linf])
				       * (wl - sixs_disc.wldis[linf]) / (sixs_disc.wldis[lsup] - sixs_disc.wldis[linf]);
			}
			else
			{
				alphaa = (float)(log(sixs_disc.roatm[2][lsup] / sixs_disc.roatm[2][linf]) / coef);
				betaa = (float)(sixs_disc.roatm[2][linf] / pow(wlinf,alphaa));
				is.roaero = (float)(betaa * pow(wl,alphaa));
			}
		}
	}

	alphar = (float)(log(sixs_disc.trayl[lsup] / sixs_disc.trayl[linf]) / coef);
    betar = (float)(sixs_disc.trayl[linf] / pow(wlinf,alphar));
    is.tray = (float)(betar * pow(wl,alphar));
    
	if (idatmp != 0)
	{
		alphar = (float)(log(sixs_disc.traypl[lsup] / sixs_disc.traypl[linf]) / coef);
        betar = (float)(sixs_disc.traypl[linf] / pow(wlinf,alphar));
        is.trayp = (float)(betar * pow(wl,alphar));
	}
    else is.trayp = 0;

	if(iaer != 0)
	{
		alphaa = (float)(log(sixs_aer.ext[lsup] * sixs_aer.ome[lsup] / (sixs_aer.ext[linf] * sixs_aer.ome[linf])) / coef);
		betaa = (float)(sixs_aer.ext[linf] * sixs_aer.ome[linf] / pow(wlinf,alphaa));
		is.tsca = (float)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
		alphaa = (float)(log(sixs_aer.ext[lsup] / sixs_aer.ext[linf]) / coef);
		betaa = (float)(sixs_aer.ext[linf] / pow(wlinf,alphaa));
		is.taerp = (float)(taer55p * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
		is.taer = (float)(taer55 * betaa * pow(wl,alphaa) / sixs_aer.ext[3]);
	}

	float drinf = sixs_disc.dtdif[0][linf] + sixs_disc.dtdir[0][linf];
	float drsup = sixs_disc.dtdif[0][lsup] + sixs_disc.dtdir[0][lsup];
	alphar = (float)(log(drsup / drinf) / coef);
	betar = (float)(drinf / pow(wlinf,alphar));
	is.dtotr = (float)(betar * pow(wl,alphar));
	float dtinf = sixs_disc.dtdif[1][linf] + sixs_disc.dtdir[1][linf];
	float dtsup = sixs_disc.dtdif[1][lsup] + sixs_disc.dtdir[1][lsup];
	alphac = (float)(log((dtsup * drinf) / (dtinf * drsup)) / coef);
	betac = (float)((dtinf / drinf) / pow(wlinf,alphac));
	float dtotc = (float)(betac * pow(wl,alphac));
	float dainf = sixs_disc.dtdif[2][linf] + sixs_disc.dtdir[2][linf];
	float dasup = sixs_disc.dtdif[2][lsup] + sixs_disc.dtdir[2][lsup];

	if(iaer != 0) 
	{
		alphaa = (float)(log(dasup / dainf) / coef);
		betaa = (float)(dainf / pow(wlinf,alphaa));
		is.dtota = (float)(betaa * pow(wl,alphaa));
	}

	is.dtott = dtotc * is.dtotr;
	float urinf = sixs_disc.utdif[0][linf] + sixs_disc.utdir[0][linf];
	float ursup = sixs_disc.utdif[0][lsup] + sixs_disc.utdir[0][lsup];
	alphar = (float)(log(ursup / urinf) / coef);
	betar = (float)(urinf / pow(wlinf,alphar));
	is.utotr = (float)(betar * pow(wl,alphar));
	float utinf = sixs_disc.utdif[1][linf] + sixs_disc.utdir[1][linf];
	float utsup = sixs_disc.utdif[1][lsup] + sixs_disc.utdir[1][lsup];
	alphac = (float)(log((utsup * urinf) / (utinf * ursup)) / coef);
	betac = (float)((utinf / urinf) / pow(wlinf,alphac));
	float utotc = (float)(betac * pow(wl,alphac));
	float uainf = sixs_disc.utdif[2][linf] + sixs_disc.utdir[2][linf];
	float uasup = sixs_disc.utdif[2][lsup] + sixs_disc.utdir[2][lsup];
	is.utott = utotc * is.utotr;

	if(iaer != 0)
	{
		alphaa = (float)(log(uasup / uainf) / coef);
		betaa = (float)(uainf / pow(wlinf,alphaa));
		is.utota = (float)(betaa * pow(wl,alphaa));
	}

	float arinf = sixs_disc.sphal[0][linf];
	float arsup = sixs_disc.sphal[0][lsup];
	alphar = (float)(log(arsup / arinf) / coef);
	betar = (float)(arinf / pow(wlinf,alphar));
	is.asray = (float)(betar * pow(wl,alphar));
	float atinf = sixs_disc.sphal[1][linf];
	float atsup = sixs_disc.sphal[1][lsup];
	alphac = (float)(log(atsup / atinf) / coef);
	betac = (float)(atinf / pow(wlinf,alphac));
	is.astot = (float)(betac * pow(wl,alphac));
	float aainf = sixs_disc.sphal[2][linf];
	float aasup = sixs_disc.sphal[2][lsup];

	if(iaer != 0)
	{
		alphaa = (float)(log(aasup / aainf) / coef);
		betaa = (float)(aainf / pow(wlinf,alphaa));
		is.asaer = (float)(betaa * pow(wl,alphaa));
	}
}
